Numerical performance of the parabolized ADM (PADM) 
formulation of General Relativity 
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In a recent paper [H the first coauthor presented a new parabohc extension (PADM) of the stan- 
dard 3+1 Arnowitt, Deser, Misner formulation of the equations of general relativity. By parabolizing 
first-order ADM in a certain way, the PADM formulation turns it into a mixed hyperbolic - second- 
order parabolic, well-posed system. The surface of constraints of PADM becomes a local attractor 
for aU solutions and aU possible well-posed gauge conditions. This paper describes a numerical im- 
plementation of PADM and studies its accuracy and stability in a series of standard numerical tests. 
Numerical properties of PADM are compared with those of standard ADM and its hyperbolic Kid- 
der, Scheel, Teukolsky (KST) extension. The PADM scheme is numerically stable, convergent and 
second-order accurate. The new formulation has better control of the constraint-violating modes 
than ADM and KST. 
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I. INTRODUCTION 

For years people have tried to obtain analytic solutions 
of the complex field equations of Einstein's general theory 
of relativity (GR) . Apart from few cases where symmetry 
is invoked, it is almost impossible to analyze the compli- 
cated dynamics in the strong gravitational field regime 
as described by GR. Approximation methods have been 
developed over the course of time, but the most promis- 
ing tool for tackling problems such as gravitational waves 
arising from binary black hole (BBH) or binary neutron 
star mergers, gravitational collapse etc., is numerical rel- 
ativity. 

In the past few years remarkable progress has been 
made towards achieving long term and stable evolution of 
the Einstein equations. Recentlyj)articular cases of the 
the BBH problem were solved [f, U, i, H, H, 0, Hi • Despite 
this remarkable achievement, the general problem of long 
term and stable evolution of the GR equations remains 
open and there is still much work left to be done. There 
is no theory or prescription to chose what formulation(s) 
and what gauge conditions are suitable for the numeri- 
cal solution of a given problem. For example, it is well 
known that the Baurngarte, Shapiro, Shibata, Nakamura 
(BSSN) formulation Q, although successful with BBHs, 
has not been successful with certain spacetimes |1Q] . In 
addition to that, there does not seem to exist a definitive 
explanation of why the approaches of 0, d, 0] perform so 
well when contrasted to previous efforts. Furthermore, 
there are astrophysical and theoretical problems of great 
interest for which the formulations above have not been 
applied yet and it is not known whether they will prove 
successful in such cases. Such problems are the study 
of the internal structure of black holes and astrophysical 
phenomena, where except for black holes matter is also 
involved. 

The numerical integration of the Einstein equations is 
not an easy task because the computations can become 



unstable and an exponential blow up of the numerical er- 
ror may occur, even when the formulation employed ad- 
mits a well-posed initial value problem. If the numerical 
techniques employed and gauge and boundary conditions 
chosen do not suffer from pathologies, perhaps the most 
important source that can potentially lead to instabilities 
during a free evolution is the growth of the constraint 
violating modes. Over the years several methods have 
been proposed to deal with this last type of instability 
[ll, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 2l,[2|. 

Of all these approaches, the one which has attracted 
most attention in recent years takes advantage of the 
fact that in the ideal case where the constraint equations 
are satisfied, one has the freedom to add combinations 
of the constraint equations to the right-hand-side (RHS) 
of the evolution equations of a given formulation. By 
virtue of this freedom, it is also possible to introduce 
in the system of evolution equations terms which act as 
"constraint drivers" and turn the constraint surface into 
an attractor. This technique is nowadays usually termed 
as "constraint damping" . 

The goal of many formulations in numerical relativ- 
ity has been to incorporate such drivers in symmetric or 
strongly hyperbolic systems without changing the prin- 
cipal part of the evolution equations [H, However, 
not all formulations with such drivers have been as suc- 
cessful as the generalized harmonic decomposition Q]. 
A possible explanation is that the constraint satisfying 
modes evolve differently with different but equivalent 
formulations. Another explanation may be that some 
formulations may require more efficient damping of the 
constraint- violating perturbations that are present in nu- 
merical simulations and lead to all sort of instabilities. 

One way to achieve efficient damping is to construct 
formulations of GR that under free evolution force 
all constraint violating modes to evolve according to 
parabolic equations. Parabolic equations are known for 
their damping and smoothing properties [23| and this 
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is of extreme importance in numerical relativity. In [l| 
this goal was achieved by the construction of an evolu- 
tion system based on the first-order form of the stan- 
dard Arnowitt, Deser, Misner (ADM) formulation ^2^, 
through the addition of appropriate combinations of the 
derivatives of the constraints and the constraints them- 
selves at the RHS of the ADM evolution equations. We 
call this evolution system the Parabolized ADM (PADM) 
formulation throughout this work. It was shown in [l| 
that the evolution of the constraint equations with the 
PADM formulation are second-order parabolic, indepen- 
dently of the gauge conditions employed. This in turn 
implies that the constraint surface becomes a local at- 
tractor. It was finally proved, that the PADM system 
satisfies the necessary conditions for well-posedness and 
based on the results of [26| an argument, which indicates 
strong evidence that the PADM system admits a well- 
posed initial value problem, was given. 

The purpose of this work is to describe a numerical 
scheme for solving the PADM equations, test the accu- 
racy and stability of the PADM system, and compare 
the PADM formulation with the first-order ADM and 
the Kidder, Scheel, Teukolsky (KST) [Hi formulations. 

The first-order ADM system in conjunction with har- 
monic or "1 -f log" slicing is strongly hyperbolic for the 
set of one-dimensional solutions considered in this work. 
Therefore, we choose the first-order ADM formulation as 
a basis for comparison. We choose the KST formulation 
because it is strongly hyperbolic and because the KST 
and the PADM systems are both extensions of the ADM 
evolution system, and we want to study how different 
extensions of the same base system perform numerically. 
Furthermore, the only method, to our knowledge, which 
has been employed to integrate the KST evolution equa- 
tions is pseudo-spectral methods. This work implements 
and demonstrates the numerical performance of the KST 
formulation with finite difference methods. 

The comparison between the three formulations is car- 
ried out in a series of standard one- and two-dimensional 
tests, usually referred to as the "Apples with Apples" 
tests or the "Mexico City" tests There are four ba- 
sic tests: a) The evolution of small initial noise on ffat 
spacetime data, b) the evolution of one-dimensional and 
two-dimensional gauge waves, c) the evolution of a small 
amplitude ID and 2D gravitational waves and d) the evo- 
lution of polarized Gowdy waves. 

This paper is organized as follows. In section [TTl we 
briefly describe the first-order ADM, KST and PADM 
formulations. In section IIIII we describe the numerical 
scheme we developed. In section lTVl we describe the stan- 
dard tests, present the results of the numerical simula- 
tions we carried out and compare the numerical perfor- 
mance the aforementioned formulations. In section |V] 
we study the convergence of the numerical schemes. We 
conclude this paper in section IVTl 



II. FORMULATIONS 

The four-dimensional theory of general relativity can 
be cast into a 3-1-1 decomposition of spacetime [23 by 
assuming that the spacetime can be foliated by a one- 
parameter family of spacelike hypersurfaces. The space- 
time metric is then written in the following form 

ds^ = -a^dt^ + 7y (da;* + l3'dt){dx^ + jS^dt) (1) 

where 7^ is the positive definite 3-metric on the t — 
const, hypersurfaces, a is the lapse function, the shift 
vector, and are the spatial coordinates, i — 1,2,3. 

A. The first-order ADM formulation 

The ADM formulation and consist of two subsets of 
equations. The first subset is that of the evolution equa- 
tions, which describe how the dynamical variables evolve 
in time. The second subset is that of the constraint equa- 
tions, which have to be satisfied for all times. The stan- 
dard second-order ADM formulation 2S] has as dynami- 
cal variables the 3-metric 7^ and the extrinsic curvature 
Kij of the 3D spacelike hypersurfaces. 

The first-order ADM formulation is derived from the 
second-order one [l| by introducing additional dynamical 
variables 

Dkij = dkjij , (2) 

and then deriving the evolution equations for Dkij ■ The 
dynamical variables of the first-order ADM system evolve 
according to the following equations 

do-fij = -2aKij, (3) 

doK,, = - V.Vj-a -I- a {R,, + KK,, - 27™"X,™ifj„) , 

(4) 

doDkij = -2adkKij - 2Kijdka, (5) 

where 7^ is the three-metric, Kij and K = j"^"'Kmn 
are the extrinsic curvature and its trace respectively, Vi 
is the covariant derivative operator associated with the 
three-metric, Rij is the Ricci tensor associated with the 
three-metric, and do = dt — £ p, with £^ the Lie deriva- 
tive along the shift vector /3'. In equation all par- 
tial derivatives of the three-metric have been replaced by 

The Lie derivatives of the dynamical variables are 

i^/37»j = V,/3j -f Vj/3„ (6) 

£pK,, = (V,DK^, + {Vjr)Krm + r^mK^,, (7) 
and 

£l3Dkij — (3"^dmDkij + DmijdkP"^ 
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The set of constraint equations is 

n= R + K^- if™„X"" = 0, (9) 

= VmK"'^-y^K = 0, z = 1, 2, 3, (10) 

Cki] = dk-fij - Dkij = 0. (11) 

where R is the trace of the three-Ricci tensor. Just hke 
in the evolution equations, all partial derivatives of the 
three- metric in ([9)) and (jlOp are replaced by Dkij- Ti. 
is the Hamiltonian constraint, are the momentum 
constraints and Ckij are new constraints due to the in- 
troduction of Dkij . The equations as presented above do 
not include matter terms, because in this work we focus 
on vacuum solutions of the Einstein equations. 

At this point we must make an important remark con- 
cerning the geometrical interpretation of the Ckij con- 
straints. The covariant derivative of any purely spatial 
tensor Aij is 

Vfev4y — dkAij — V^ikAsj ~ T^kjA-is, (12) 

where, in the context of the first-order ADM formulation, 
T'^ ij is given by 

r^-^7'"(i?.» + A,n~i?™,). (13) 

The covariant derivative of the three-metric must be zero, 
but by virtue of and (fT^ it is straightforward to 
show that 

Vfc7jj = Ckij- (14) 

Therefore, the geometrical interpretation of the Ckij 
constraints is that the covariant derivative is metric- 
compatible, if and only if Ckij = 0. 

To close the system, equations - (0) must be supple- 
mented with the gauge equations for the lapse function 
and the shift vector. Following [27i] . we write them in a 
general form as 

Fa(^^,a,l3\dba,dbdca, ...,dt(3\ ■■■jij,dt-fij, ...^ = 0, 
a,5,c= 0,...,3, i,j = l,2,3. 

(15) 

In this work we limit our consideration to algebraic 
gauges where the lapse function is of the form 

a — 0(7), (16) 

where 7 is the determinant of the three- metric jij, and 
the shift vector /3* may be either constant or a fixed func- 
tion of the spacetime coordinates {t, x^). In [2^,[23| it was 
shown that makes the evolution well-posed on the 
surface of constraints, \i A = dhia/ dXn^ > 0. Working 
with (fTO|) . we have the flexibility to use either "l-flog" 
slicing 

a = l + ln7 (17) 



or a densitized lapse 

a = Q7", (18) 

where Q is constant or a fixed function of t, and a the 
densitization parameter. If p7|l is used, A ^ \/a. If 
PS)) is employed, A = a. 

B. The KST formulation 

The strongly hyperbolic, four-parameter KST modifi- 
cation to the first-order ADM formulation |28] is 

dtK.j = (ADM) + pa-i.^n + i^a7"''C„(y)fc, (19) 

^tDk^J - [ADM) + f]a-ik(:Mj) + xal^]Mk, (20) 

where [ADM] stands for the RHS of the first-order ADM 
evolution equations and p, tp^ rj and x £^re four parame- 
ters of the formulation. In addition to those modifica- 
tions, Kidder, Schccl and Teukolsky used a densitized 
lapse gauge (fT8)) . 

In this work we use the KST system in conjunction 
with the more general gauge condition p6p . We do so 
because in fl*] it was shown that the KST formulation 
with (|16|) is strongly hyperbolic, if the KST parameters 
are 

V=l, ^ = '1^ X = -| Pt^O. (21) 

It was also shown that for the particular choice of p = 
— 1/3 the number of flat space modes which violate the 
Hamiltonian constraint is smaller than for any other 
value of p. Thus, the values of the KST parameters that 
we will use throughout this work are those in (PT|) with 
P = -l/3. 

C. The PADM formulation 

The PADM system is obtained from the first-order 
ADM formulation by addition of constraints and their 
derivatives to the RHS of the ADM evolution equations. 
It has six parameters and is given by 

9t7y = (ADM) + Xr''^bCa^J, (22) 
dtK.j = (ADM) + Hvl^'daMb + ed^.Mj) , (23) 

^tDk^■J = {ADM) + er''^aCbk^o + ^,l^o^k^i + C^Ck^, 

(24) 

where A, 0, 0, e, ^, C are the six parameters of the formula- 
tion. We refer to the added terms of the PADM formula- 
tion as the constraint driver. It was shown in [l| that for 
gauges which do not introduce second-order derivatives 
of the dynamical variables in the evolution equations the 
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PADM system satisfies the Petrovskii condition for well- 
posedness [2^, provided that 

A > 0,e > 0,^ < 0,(/)< 0,6* > 0. (25) 

Therefore, in this work we use the PADM equations in 
conjunction with the algebraic gauge condition (|16p . 

When the aforementioned conditions are satisfied, the 
PADM system has the following properties: a) The 
PADM evolution equations can be classified as a set of 
mixed hyperbolic - second-order parabolic quasi-linear 
partial differential equations. The parabolic character of 
the equations can be most easily seen by the evolution 
equation of the 3-metric, whose principal part is given 
by 

-xr'dadbj^j, (26) 

where ~ implies equal to the principal part. 

b) The evolution equations of the constraint variables 
become second-order parabolic partial differential equa- 
tions (PDEs) independently of the spacetime geometry 
and the gauge conditions employed. Therefore, the con- 
straint propagation equations admit a well-posed Cauchy 
problem themselves, c) Because of the parabolic struc- 
ture of the evolution equations of the constraints, the 
constraint surface becomes a local attractor. All small 
amplitude, high-frequency constraint-violating perturba- 
tions are exponentially damped in time as exp(— Ak^^), 
where k is the magnitude of the wavevector of the per- 
turbations. This in turn implies that the hazardous high 
frequency perturbations must be damped very efficiently. 

In this work we use the PADM formulation with the 
following set of parameters 

( = 1,A = e = 61 = -2^ = -20 = 0.02 (27) 

in all our simulations. We do this for two reasons: a) It 
was shown in [IJ that for perturbations about flat space, 
the Fourier transformed operator of the evolution equa- 
tions possesses a complete set of eigenvectors if 

A = e = -2C = -2<j) ^ 9/2 (28) 

and b) when using explicit numerical schemes to carry 
out the integration of the PADM evolution equations, the 
values of the PADM parameters given in (^7)) are small 
enough to allow for large enough time-steps, while the 
damping properties of the formulation are still present. 

Finally, we note that care must be taken in the choice 
of the parameters of the PADM system for backwards in 
time evolutions, e.g, the collapsing Gowdy spacetimes. 
The backwards in time Cauchy problem for a parabolic 
equation is ill-posed, because of the existence of exponen- 
tially growing modes. The same is true for the PADM 
formulation if the parameters of the formulation satisfy 
(I25p . To overcome this problem we simply have to reverse 
the signs of the PADM parameters in . 



III. NUMERICAL SCHEME 

We use a finite difference method to carry out com- 
putations. For all formulations we use a staggered spa- 
tial mesh in which jij and Kij are located at the cell 
centers and D^ij are staggered half a grid point in all 
spatial directions. Moreover, all variables are defined at 
the same time layer. For ADM and KST we use a third- 
order Runge-Kutta method to do the integration in time, 
as this scheme has shown to posses desirable dispersion 
and dissipation properties compared to other commonly 
used numerical schemes [s^. We use second-order ac- 
curate centered derivative operators to calculate spatial 
derivatives and a 3rd order accurate parabolic interpola- 
tion operator whenever staggered values of the dynamic 
variables are needed. 

The PADM formulation has both hyperbolic and 
parabolic terms in the RHS of its evolution equations. 
Since we have implemented an explicit algorithm, if we 
were to deal with the integration of the PADM evolution 
equations in the same way as for the KST and ADM sys- 
tems, we could soon face strict limitations on the speed 
of the computations because of limitations on the maxi- 
mum allowed time-step. We deal with this limitation by 
combining the methods of operator split and fractional 
steps. 

The basic idea of operator splitting was originally in- 
troduced in [3l| . Here we omit all the details and present 
the technique in its simplest form. Let 



be a set of PDEs, where u is the column vector of the un- 
known variables and £ is a differential operator. Assume 
further that C can be written as a sum of two operators, 
i.e., 

Cu = Ciu + C2U. (30) 

Let us now consider the individual equations 

dtu = Ciu, dtu = C2U (31) 

and let C(Ai) be a finite difference operator for equa- 
tion (I29|) . i.e., the solution at time-step n -|- 1 is given by 
it"+i = C(Ai)u". Also, let Ci and C2 be corresponding 
finite difference operators for the individual equations of 
([3T|l . If the operators Ci and C2 are second-order accu- 
rate and stable for time-step At, then it can be shown 
that the approximation 

C(2At) = Ci(Ai)C2(Ai)C2(Ai)Ci(Ai) (32) 

provides a second-order accurate and stable scheme for 
(pg)) for a time-step 2At. 

In the discussion above we made no explicit mention 
to the exact form of the finite difference operators Ci 
and C2. Therefore, one is free to use any finite difference 
operators, as long as those are stable and second-order 
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accurate for a time-step At. This means that for any of 
the operators C i an d C2 we can use the scheme known as 
fractional steps [S^ . In this scheme instead of advancing 
the solution forward a whole time-step At, one advances 
successively q times, each with a time-step of At/ q. Thus, 
if we choose to do fractional stepping for the C2 operator, 
we can approximate C2(Ai) as 

C2 (At) = C2{At/q)C2{At/q)...C2{At/q) (33) 
g-times 

One may naturally think of the PADM formulation as 
consisting of two operators which act additively on the 
dynamical variables, such as those described by equations 
([59)1 and dSO]). The first operator is that of the RHS of the 
evolution equations of the ADM formulation and the sec- 
ond one is the constraint driver. There is no unique way 
to split the entire operator of the RHS of the evolution 
equations of the PADM formulation, but because of its 
mathematical structure and for computational efficiency, 
the most straightforward way to proceed with splitting 
is to set the £1 operator equal to the ADM operator in- 
cluding the low order terms of the constraint driver, and 
set £2 equal to the higher order terms of the constraint 
driver. We call this Li operator the hyperbolic operator 
and £2 the parabolic operator. If the splitting is done 
in this way, then the £2 operator has smaller number of 
operations on the dynamical variables compared to any 
other choice of splitting. This feature is important for 
computational efficiency as it will be evident shortly. 

For the PADM formulation, we use the techniques de- 
scribed above in the following way. We operator split 
equations as is described in the previous para- 

graph. For the time integration of the hyperbolic terms 
we use the same technique as that for the ADM and KST 
formulations. For the time integration of the parabolic 
terms we use a second-order accurate iterative Crank- 
Nicolson scheme with two iterations. Furthermore, when- 
ever it is required by the numerical stability criterion, we 
use the fractional steps method to integrate the parabolic 
terms. We do this by finding the smallest positive integer 
number p, defined by 



At„ 



At 

P 



(34) 



for which the parabolic time-step Atp satisfies the stabil- 
ity condition (see discussion on numerical stability below) 
and then proceed with the integration of the parabolic 
terms as is dictated by p3p . In equation (f34|) . At is 
the time-step of the hyperbolic terms. The advantage of 
this algorithm is that even though we may have to use 
a large p for high resolutions, and thus have to calculate 
the parabolic terms p times for each hyperbolic time-step, 
it reduces the computational cost significantly, because 
the hyperbolic terms have larger computational overhead 
than the overhead of the parabolic terms. 

An important issue which has to be addressed at this 
point is the numerical stability criterion of the algorithm 



for the PADM formulation. Since there is a hyperbolic 
and a parabolic part in the equations, there arise two 
stability conditions: a) A condition from the hyperbolic 
part, and b) a condition from the parabolic part. There- 
fore, we must always make sure that we satisfy both the 
hyperbolic and the parabolic stability criteria [3^ . 

The stability criterion for general hyperbolic equations 
with anisotropic speed of propagation of information is 



min{ 



CxAt 



CyAt 



.At 



) = cfl<ci, 



(35) 



where Cx, Cy, and Cz are the speeds of propagation of in- 
formation in the x—, y—, and z— directions respectively, 
cfl is the Courant-Friedrichs-Levy number, and ci is a 
constant which depends on the numerical scheme. For 
example, for a scalar wave equation with speed unity in 
conjunction with a 3rd-order Runge-Kutta time integra- 
tor and second-order accurate centere d sp atial derivative 
operators the constant is ci — ^3/4 |3(J |. 

For most applications in numerical relativity the speed 
of propagation of information is isotropic and the choice 



At 

— = 0.25 
Ax 



(36) 



satisfies condition (|35p . For this reason in this work we 
set the time-step for the hyperbolic part of PADM via 
((36)) . We do the same for all computations with the ADM 
and KST formulations. 

Now we turn our attention to the stability criterion of 
the parabolic part. A Von Neumann analysis of equation 
(j26p shows that in order to preserve numerical stability 
the explicit numerical scheme we use must satisfy 



AAt, 



7 



11 



-I- 



r 



22 



+ 



7 



33 



{Axy (Ay)2 (Az) 



< 



1 



(37) 



Although condition ([37]) results from the stability analy- 
sis of the principal part of the evolution equation of the 
3-metric, all our numerical experiments confirm that if 
we choose the PADM parameters according to 



A = e 



-2^ = -2(j> 



(38) 



then if ([37]) is satisfied the computations remain stable, 
and if it is violated the computations soon become un- 
stable. 

Taking all these facts into consideration, we conclude 
that if condition ([55]) is satisfied, the overall numerical 
stability criterion of the numerical scheme for the PADM 
system is 



p> 



A(7"+7''+7'') 



2Ax 



(39) 



with p = 1 the minimum value p can obtain. In the 
derivation of the last equation we combined (|M)) . (I36p . 
and §7}. 
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IV. NUMERICAL TESTS 

In this section we briefly describe the standard tests 
we used and the results of the numerical comparison of 
the three formulations presented in section|TTl As is usual 
in the literature, we set up tests which probe both the 
linear and the non-linear regime of the GR equations. 
We closely follow 10] in the setup of the tests, but also 
modify some of them in order to make the distinction 
among the three-formulations more apparent. 

For the evaluation of errors we follow [33] and check the 
accuracy of the integration with a single number given by 
the following L2 norm 



where Vol = / ^(Px, and 



(40) 



(41) 



The squared quantities are equal to the norms of the 
quantities with respect to the numerical three-metric ten- 
sor, e.g., [CkijY = C'^^Ckij, where all indices are raised 
with the numerical metric. It is clear from equation (|40|) 
that 11C]]2 — 0, if and only if all the constraints are sat- 
isfied. 11C))2 is called the "constraint energy." 

Similarly, we monitor the accuracy of the numerical 
solution as compared with the analytic one, by using the 
"error energy" which is given by 



where 



SU 



(42) 



(43) 



and where all "delta" dynamical variables constitute the 
difference between the numerical and the analytic solu- 
tion, e.g. = -^analytic _ ^numerical_ ^s in the con- 
straint energy, the indices in the error energy are raised 
with the numerical metric. The error energy vanishes, if 
and only if the numerical solution matches the analytical 
one. 

Just like in ^] , for the Gowdy wave tests we plot the 
normalized version of those quantities, i.e., ] J^Z/^J I2/I JZ^J I2 
and 11C112/119Z^112, where 



- ^J{l^,? + {K^,? + {Dk^,)\ (44) 
dU = ^{^kl^,Y + {dkK,,Y + {dsD^jY, (45) 



and 



and where the squared quantities here are with respect to 
the analytic three-metric tensor. The reason for plotting 



the normalized errors is that the error depends on the 
magnitude of the solution. 

We set up all tests on a 3D mesh with periodic bound- 
ary conditions. We perform all one-dimensional tests 
on a cell-centered domain, such that x £ [Ojl]i with 
Xn — {n — ^)Ax, where n — 1,2,... 2^ and Ax = 
Ay = Az = 2"'', where r € N. In addition, we fol- 
low [3^ and set up the domain with four points in the 
y— and z— directions. In the case of two-dimensional 
tests the domain is similar. In particular, x £ [0,1], 
y G [0, 1], Xn = {n- ^)Ax, y„ = (n- ^)Ay, where again 
?i = 1, 2, . . . 2^, Ax = Ay = Az = 2^^, and have 4 points 
in the z-direction. 



A. Robust stability test 

The robust stability test is based on small random per- 
turbations about Minkowski spacetime. The amplitude 
of the perturbations is sufficiently small so that the evo- 
lution remains in the linearized regime unless any insta- 
bilities arise. The initial data are given by 



7y 



D 



kij 



(46) 



where the random numbers {^Ij , f^j , ^kij} are uniformly 
distributed in [-lO^^Vp, lO^^Vp], where p = 2'-^ 
Note that instead of rescaling by p^, as is suggested in 
p^ . here we rescale by p. We do this because the con- 
straints for the formulations we are considering are first- 
order, and in order to achieve the same initial violation 
on the constraints we have to rescale by p. We call the 
above test the "standard" robust test. 

Another version of this test requires perturbation by 
noise not only at the initial time-step, but at every time- 
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FIG. 1: Constraint energies as functions of coordinate time t 
for the standard robust test with a = 7^''^. The letters a, b, c 
correspond to resolutions r = 4, 5, 6 respectively. Curves la, 
lb and Ic correspond to the ADM formulation, curves 2a, 2b 
and 2c correspond to the KST formulation and curves 3a, 3& 
and 3c correspond to the PADM formulation. 
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step throughout the entire simulation. We call this test 
the "strong" robust stability test. We run both versions 
of the robust test with r = 4, 5, 6 for 1000 light crossings 
or until the code crashes. 



1. Standard Robust Stability Test 

The results of this test ran with harmonic slicing are 
summarized in Figure [U where we have plotted the con- 
straint energy for the ADM, KST and PADM formula- 
tions as a function of time. The ill-posed nature of the 
ADM formulation is clearly seen in this figure . The 
simulations crash fairly quickly and most importantly at 
earlier times for higher resolutions, as expected for ill- 
posed systems. 

Figure [T] also shows that the three error energies of the 
KST runs practically overlap. This is a result of rescal- 
ing the amplitude of the random noise with resolution, 
so that the initial violation of the constraints is of the 
same magnitude for all resolutions. The KST simula- 
tions exhibit noticeable growth in the constraint energy. 
Close investigation shows that Ckij grow linearly with 
time, whereas 7i and M.i stay roughly constant through- 
out the entire simulation. Therefore, the growth of the 
constraint energy comes from Ckij ■ 

The linear growth with time of Ckij has a very simple 
explanation. The evolution equations of the Ckij con- 
straints with the KST formulation are 

doCkij = -mik(iMj) - xajijMk- (47) 

This last equation implies that any small violation of 
the Momentum constraints will be spilled into the Ckij 
constraints and cause them to evolve with time. Our 
results of the robust test with the KST formulation show 
that the RHS of equation (IT71) is roughly equal to a non- 
zero constant. This implies that the Ckij constraints must 
grow linearly with time. 

We will now show that the linear growth with time 
of Cm with the KST formulation is due to violations of 
Aix- We will do this by using our results of the r — 6 
runs of the robust test. Via a least squares fit of the 
infinity norm of Cm as a function of time we find that its 
observed growth rate is dt\Ciii\°^'"''''""^ « 6.5773 • 10~^°. 

To show that the linear growth of |Ciii|oo is due to vio- 
lations oi Aix we need to demonstrate that the observed 
growth rate agrees with the predicted growth rate from 
equation (|47p . Keeping in mind that the lapse function 
and 7ii are practically equal to unity throughout the en- 
tire run, and using ((2T|) and (|47|) we find that the absolute 
value of the time derivative of Cm is 

dt\Cin\ - l\Mx\. (48) 
5 

Using the data from the same simulation we find that 
the infinity norms of |A^2;|oo stays roughly constant dur- 
ing the entire run and that its time-average is such that 



dt\Ciii\P^'"^'''*'"' « 6.162 • 10-1°. The growth rate pre- 
dicted by equation (|T7l) agrees with the observed one 
within 6.3%. 

A consequence of this pathology of the KST system 
is that given sufficient time or stronger initial pertur- 
bations, the linear growth will eventually terminate any 
KST simulation. 
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FIG. 2: Cumulative amplitude Ak (see equation (|50p ) of the 
Fourier transform of the constraint variables as a function of 
coordinate time t, for harmonics k = 1,2,4,8,16,32. The 
plot corresponds to the run of the standard robust test using 
a = 7^''^, with the PADM formulation, for resolution r = 6. 

Figure [T] shows that the PADM formulation has dras- 
tically different behavior than the ADM and KST sys- 
tems. The constraint energy decreases very rapidly by 
many orders of magnitude and then reaches a plateau. 
This formulation damps any small violations of the con- 
straint equations and pushes the evolution back onto the 
constraint surface. The plateau can be explained as fol- 
lows. Because of roundoff errors, random noise of the 
order of 10" is introduced into the system at every 
time-step. This roundoff error noise causes violations of 
the constraints of the order of 10~^^/dx because all our 
constraints contain first-order derivatives of the dynami- 
cal variables. However, dx is roughly of the order of 10"^ 
and hence the average magnitude of the constraints when 
they reach the plateau is expected to lie at roughly 10^^*. 
This is in accordance with what we observe if Figure [1] 

The plateau occurs at the point where the damping 
of the constraints essentially balances the constraint vi- 
olations caused by the roundoff errors. There is a slight 
difference in the plateau at which the constraint energy 
fiattens out for different resolutions, ft is seen that for 
the highest resolution r — 6 (curve 3c), the plateau lies 
higher than those for resolutions r — A and r — 5 (curves 
3b and 3a respectively). This is because the roundoff 
errors are on average of the same magnitude and hence 
the roundoff error perturbations are not scaled as 1/p, as 
is the case with the initially added random noise. Con- 
sequently, the plateau of simulations of the robust test 
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(a) 

FIG. 3: Constraint energies as functions of coordinate time t 
for the standard robust test in l+Zog slicing. The letters a, 6, c 
correspond to resolutions r = 4, 5, 6 respectively. Curves la, 
16 and Ic correspond to the ADM formulation, curves 2a, 26 
and 2c correspond to the KST formulation and curves 3a, 36 
and 3c correspond to the PADM formulation. 

with the PADM system wiU he higher with increasing 
resolutions. 

In Figure [2] we explore the damping properties of the 
PADM formulation. The plot in the figure was created 
as follows. Since the simulation is quasi one-dimensional, 
we focus on a one-dimensional line located along x — 
[0,1],?/ = /S.y/2^z = Az/2. Let the constraint variables 
be denoted by Ci, where i = 1, 2, . . . , 18. Along this line 
and at a given time t we perform a discrete Fourier trans- 
form of the N = 2"^ values of all the constraint variables 
Ci,n, n — 1,2, . . . , N , so that 

N 

Aa- = ^C,,„e-'^'=", k = Q,l...,N-l (49) 

where Ai^k denotes the fourier harmonic k of the con- 
straint Ci. For a given harmonic, we then sum the am- 
plitudes of all the Fourier pairs of the constraints, i.e., 

18 

Ak=Y,\^^M- (50) 

i=l 

Ah is zero, only if all constraints have zero contribution 
to the specific harmonic k. 

The quantity Ak as a function of time is what is plotted 
in Figure [2] for the standard robust test with the PADM 
formulation. Ak provides us with a measure which de- 
scribes the cumulative amplitude of a specific Fourier 
mode (fc) over all the constraints. The set of A}, also car- 
ries the information of whether any instability arises and 
on which end of the frequency spectrum it occurs. Figure 
[5] shows that the PADM formulation damps the highest 
frequencies (higher k) the most. This is precisely what 



was predicted in [l| for the PADM system. It is worth 
noting that the Nyquist frequency (32nd harmonic) con- 
straint violations are damped away in less than one light 
crossing time and all the constraint violating modes are 
damped away after approximately 25 light crossings. 

Finally, we also carried out the standard robust stabil- 
ity test with 1 -|- log slicing and with a ct = 1.0 densitized 
lapse and we find similar behavior for all formulations. 
For brevity we only show the results for 1 -I- log slicing 
(Figure [3]). It is again worth noting that the PADM 
formulation damps the constraints, indicating that the 
damping properties of the formulation are independent 
of the gauge choicej_as long as the gauge chosen is well- 
posed according to [2^). This again is in excellent accor- 
dance with the predictions of 



2. Strong Robust Stability Test 

This test is much more challenging than the standard 
robust test and is much closer to reality because in nu- 
merical simulations, truncation error will unavoidably be 
introduced into the solution of the PDE system at every 
time-step. 

Simulations of this test with the ADM formulation 
crash almost immediately. For this reason we only show 
the results for the KST and PADM formulations. In Fig- 
ure [4] we show the constraint energies resulting from the 
strong robust tests for the KST and PADM systems for 
various noise amplitudes. All plots in this figure corre- 
spond to resolution r — 6. For the KST formulation 
we see that at amplitude A = 10~^° (curve la), the 
constraint energy grows more rapidly than in the stan- 
dard robust test (Figure [1]). Nevertheless, the simulation 
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FIG. 4: Constraint energies as functions of coordinate time t 
for the strong robust test, with a = 7^''^ and resolution r = 6. 
The letters a,b,c correspond to noise amplitudes A — 10^^'', 
A = 10~* and A = 10~^ respectively. Curves la and 16 
correspond to the KST formulation. Curves 2a, 26 and 2c 
correspond to the PADM formulation. 
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FIG. 5: Errors as functions of coordinate time t for the ID gauge wave, with A — 0.1 and resolution r — 7. On the left panel 
(a) we show the constraint energies and on the right panel (b) we show the error energies for the ADM, KST and PADM 
formulations. The PADM formulation damps the constraints within three light crossing times and reaches the plateau value 
||C||2 = 3.11 ■ 10"**, while the ADM constraint energy at t = 1000 is equal to ||C||2 = 4.15 ■ lO"'^. 



can be run for 1000 light crossing times. However, it 
is clear that given enough time, the constraint energy 
will eventually grow so large that the simulation will be 
terminated. This can be seen more clearly in curve lb, 
which shows the constraint energy for noise amplitude 
A = 10~^. It is evident that in this case the KST runs 
manage to complete only approximately 600 light cross- 
ings before the simulation is so polluted that it crashes. 

The results of the PADM runs are dramatically differ- 
ent. Figured] shows no sign of increase of the PADM con- 
straint energy for the first 1000 light crossings. Naturally 
we do not see the damping behavior to the level shown 
in Figure [TJ because strong noise is being pumped into 
the system at every time-step. However, it is impressive 
that even with random noise of amplitude A = 10~^ the 
constraint energy does not grow throughout 1000 light 
crossings and the simulations show no signs of patholo- 
gies. 



B. Gauge wave test 



The metric of the one-dimensional gauge wave test jlC 



IS 



where 



ds^ = -hdr + hdx^ + dy^ + dz\ (51) 



/i = 1 + yl sin 



2tt{x - t) 



(52) 



where A is the amplitude of the gauge wave and d the 
size of the evolution domain. For the first-order ADM 



formulation (|5ip . ([5^ is equivalent to 

'2Tr{x - t) 



7ii =1 + A sin 



All = ; cos 



d-fi 



d 

2-k(x - t) 



(53) 



722 = 733 = 1, and all other dynamical variables are zero. 
A further coordinate transformation 



2nA /2Tr(x-t) 
Dm = — — cos ' 



1 



V2 



{x'-y'), y 



1 



V2 



(54) 



yields the two-dimensional diagonally traveling gauge 
wave. 

The gauge wave stability test has proven to be one 
of the most challenging tests for most evolution schemes 
mainly due to the existence of a one-parameter family of 
exponential, harmonic gauge solutions [35| which corre- 
spond to the metric 



ds' 



e^*/i(-dt2 + dx^) + dy^ + dz^ 



(55) 



for arbitrary ^. The gauge wave test ((5T|) corresponds to 
fi ~ 0, but in any numerical implementation of harmonic 
gauge conditions the exponential modes, /i 7^ 0, are likely 
to be excited by truncation error. This eventually causes 
the numerical solution to veer off the analytic one. 

The growth of these exponential gauge modes can be 
avoided by the incorporation of a set of semi-discrete 
conservation laws into the principal part of the evolu- 
tion equations of a harmonic formulation (ssj . However, 
these conservation laws are specific to the gauge- wave so- 
lution and it is not known whether similar conservation 
laws exist for arbitrary spacetimes. 
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The lapse function for the gauge wave test is harmonic 
and we use a — 7^/^ in our simulations. Since these 
coordinates are harmonic and we do not incorporate any 
conservation laws in our scheme, we do expect to excite 
the exponential gauge modes. We run both the ID and 
the 2D gauge wave tests with amplitudes A = 0.01,0.1 
and resolutions r — 5,6,7. The ID gauge wave is run 
for 1000 light crossings or until the code crashes and the 
2D gauge wave is run for 100 light crossings or until the 
code crashes. Below we describe the results of this test 
for each formulation separately. 



1. The low amplitude ID gauge wave 

The A = 0.01 ID gauge wave presented no appar- 
ent difficulty for any of the formulations used in this 
work. For brevity, we only describe the results of this 
test because the high amplitude gauge wave, discussed 
directly below, is much more interesting. The error en- 
ergies as functions of time, with all three formulations, 
completely overlap. However, while the KST and ADM 
constraint energies overlap, the PADM constraint energy 
lies at least two orders of magnitude below those two. 
The PADM formulation damps the constraints very fast 
and reaches the plateau value ||C||2 = 3.17 • 10^^ within 
three light crossing times. 



2. The high amplitude ID gauge wave 

Our results for the high amplitude ID gauge wave are 
summarized in Figure [5l where we show the error and 
constraint energies as functions of the coordinate time t 
for resolution r = 7. 

a. The behavior of the ADM formulation: Figure 



5(a) shows that the constraint energy in the high am- 
plitude gauge wave runs with the ADM formulation does 
not grow. However, Figure [5 (b) | shows that the ADM er- 
ror energy grows quickly with time and overlaps with the 
PADM error energy. For both systems, the source of this 
growth is the excitation of exponential harmonic gauge 
modes. We do not study here this subject because we 
address it thoroughly below in the context of the PADM 
results of the ID gauge wave test. 

b. The behavior of the KST formulation: Figure [5] 
show that the KST system cannot pass the high ampli- 
tude gauge wave test because the computations do not 
complete 1000 light crossing times. The KST constraint 
energy suddenly begins to grow exponentially at i « 160, 
and later on at t « 200 the KST error energy blows up. 

We will now show that the unexpected behavior of the 
KST system in the simulations of the A — OA gauge wave 
is due to excitation of unwanted degrees of freedom by vi- 
olations of the momentum constraints. If there is no noise 
introduced in an evolution system, the only evolving dy- 
namical variables for the one dimensional gauge wave 
must be 711, -R'n, and Dm. This implies that the only 
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FIG. 6: Infinity norms of the Cm , H, A4x and C122 constraints 
as functions of the coordinate time t. The plot corresponds 
to the KST run for the gauge wave test, with A = 0.1 and 
resolution r — 7. 



constraint that evolves is Cm. Furthermore, Ti, and Aix 
must evolve due to roundoff error, but must be of the or- 
der of 10~^^. However, this is not the case with the KST 
system. It is straightforward to show by use of (f20|) that 
the KST equations excite D122, -D133, D212, and D313 be- 
cause Aix is not exactly satisfied. These dynamical vari- 
ables in turn excite 722,733,712,713, -^^22, i^33, -f^i2, ^^^13 
and eventually all the constraints these variables are in- 
volved in, i.e., C122, C133, C112, C113, Mi, and H, grow 
rapidly with time. The constraint violations soon become 
very large that a non-linear instability develops and the 
evolution is terminated. 

The excitation of unwanted degrees of freedom can 
be seen in Figure [SI where we show the infinity norms 
of Cm, C122, Mx and Ti. as functions of the coordi- 
nate time t. An instability causes C122, Mx and Ti. to 
grow exponentially in time from very early on in the evo- 
lution. When the violations of C122, Mx and Ti grow 
large enough, the solution gets spoiled and at t ~ 200, 
the numerical solution is so polluted that the simulation 
crashes. 

To prove that that the instability of the evolutions of 
the gauge wave test with the KST formulation is caused 
by the excitation of -D122, -0133, D212 and -D313 by round- 
off error in Mx, we also ran a "partially constrained" 
evolution of the gauge wave with the KST system, where 
we set the D122, -D133, D212 and D313 equal to the exact 
values these should have at every time-step, i.e., 0. The 
result of these runs was that the KST formulation can 
evolve the gauge wave for 1000 light crossing times with- 
out any problems. Furthermore, we find that the KST er- 
ror and constraint energies completely overlap with those 
of the ADM formulation. We do not show plots of these 
simulations, but any plots can be made available by the 
authors upon request. 

These results indicate that simulations of the gauge 
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wave with formulations which add multiples of the con- 
straints to the RHS of evolution equations are likely to 
be terminated by the same type of instability that termi- 
nates corresponding simulations with the KST formula- 
tion. An example of such formulation is the BSSN system 
and it is known that until now no BSSN based numeri- 
cal scheme has been able to giv e satisfactory gauge wave 
simulations with A = 0.1 |3^. [sgI . [37] , Moreover, the 
hyperbolic system studied in [38| demonstrated similar 
behavior. We do not analyze this subject any further, as 
it is beyond the scope of this paper and will be addressed 
in a future work. 

c. The behavior of the PADM formulation: Figure 



5(a) shows that the PADM constraint energy does not 
grow with time and that it is over an order of magni- 
tude less than the corresponding ADM constraint energy. 
However, Figure 5(b) shows that the PADM error energy 



grows quickly in time. We will now show that the source 
of this growth in the gauge wave simulations is the exci- 
tation of exponential harmonic gauge modes ([55]) . 

To obtain direct evidence that exponential harmonic 
gauge modes are excited during a simulation we must 
show that at any given time t 



Ixx = f{t)lll, 



(56) 



where ^xx is the numerical solution corresponding to 711 
of equation (j53p . and where f{t) is a normalization fac- 
tor to be determined. We do not assume that f{t) — e^* 
with fixed /i because truncation errors are continuously 
introduced into the numerical solution and hence differ- 
ent values of /i are likely to be excited as the evolution 
proceeds. 

We will now show how one can calculate f{t) nu- 
merically. The spatial average of 711 is unity, because 
sin(27r(x — t)/d)dx = 0. Therefore, for an exponen- 
tial harmonic gauge mode, ^ii{^,,x,t) ~ e''*7ii(a;, i), the 
analytical counterpart fa{t) oi f{t) is 



fa{t) 



-fiiip., x,t)dx. 



(57) 



The result of the integration of equation ([57)1 is fa {t) = 
e'^* as expected. To find f{t) all we must do is to replace 
7ii(/x,a;,t) in (|57p by ^xx- fit) must be different for dif- 
ferent resolutions, because ^xx changes with resolution. 
For this reason we denote the normalization factor of 
at the n— th time-step by /" and we calculate it as 



1 



(58) 



where 7"^ , denotes the value of the numerical solution 
at the i-th grid point and time tn- 

To demonstrate that equation ([5S|) holds true in our 
simulations of the gauge wave, we use the data of the 
gauge wave runs with A = 0.5. In Figure [7] we show a 
snapshot ofjxx, Ixx/fr ^^'^ 7ii at i = 300. The continu- 
ous curve corresponds to the numerical solution with the 



PADM formulation and the dashed curve to the analyti- 
cal one. The rescaled numerical solution corresponds to 
the dotted curve, but is essentially overlapping with the 
exact solution. We find that the rescaled numerical solu- 
tion overlaps with the analytic one not only ai t — 300, 
but also at all times. This proves that ([56]) holds true in 
our simulations and in turn provides us with direct evi- 
dence that the main source of the growth of the error is 
the excitation of exponential harmonic gauge modes by 
truncation error. The same analysis as presented above 
shows that the same source of error causes the growth of 
the PADM error energy seen in Figure |5(b)[ 

We now turn our attention to the properties of function 
f{t) and we will show that in our simulations f{t) = 
e°'* , a > 0. Figure [5] shows the natural logarithm of /" 
as a function of time for the gauge wave runs with A = 0.5 
and resolutions r = 7, 8. If truncation errors excited a 
given of the solution ([55]) , the curve i — In /" should be 
a straight line with the slope equal to ^. However, the 
growing slope of the curves in Figure [5] indicates that the 
normalization factor grows faster than exponentially. If 
we interpret the local in time slope of these curves as the 
value of /I, the former implies that /i grows as time goes 
by 

In order to find how fi evolves with time we assume 
that f{t) = e'^* with time- varying /i. If this is the case, 
fi — hif(t)/t. Thus, to find ^{t) we need to study how 
In f^/t evolves with time. We do not show t — In /" /t 
plots, but our results indicate that ln/"/t grows lin- 
early with time. This implies that n(t) — a ■ t and thus 
/(<) = e°* . For example, by using the r = 8 data we 
find via a least squares fit that a = 1.05915 • 10~^. As 
the resolution increases, a decreases. This is also clear 
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FIG. 7: Snapshot of the solution of the ^xx component for the 
ID gauge wave, with A — 0.5, resolution r = 8 and t — 300. 
The continuous curve corresponds to the PADM numerical 
solution, while the dashed curve is the exact solution. The 
dotted curve corresponds to the normalized PADM numeri- 
cal solution, which is essentially overlapping with the exact 
solution. 
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FIG. 8: Normalization factor /" as a function of time. The 
plot corresponds to the PADM simulations of the ID gauge 
wave, with amplitude A = 0.5 and resolutions r = 7, 8. 



in Figure [H because at any given time < /" . This 
is consistent with the fact that with increasing resolu- 
tion the truncation errors become smaller and hence the 
solution approaches fi = 0. 

The existence of exponential harmonic gauge modes is 
known to be a property of harmonic coordinates only. 
Thus, we anticipate that simulations of the gauge wave 
would not demonstrate the same error growth in other 
gauges. This analysis will the subject of future work 
and for now we conclude that the performance of the 
PADM formulation in the ID gauge wave test is better 
than those of the ADM and KST formulations, since the 
overall PADM error is less than the overall ADM error 
and much less than overall KST error. 



The low amplitude 2D gauge wave 



The simulations of the 2D gauge wave with A — 0.01 
are not as interesting as those of the 2D strong gauge 
wave. For this reason we only describe the results of the 
weak 2D gauge wave here. The PADM and KST for- 
mulations have no problem completing 100 light crossing 
times without any pathologies in the constraint error. 
The constraint energies of the simulations with PADM 
are damped and lie lower than the simulations with KST. 
On the other hand, both the PADM and KST systems 
experience growth in the error energy due to the excita- 
tion of exponential harmonic gauge modes. In contrast to 
the KST and PADM results, the ADM formulation can- 
not evolve the initial data for more than 66 light crossing 
times. 



4- The high amplitude 2D gauge wave 

Our results of the strong two-dimensional gauge wave 
are summarized in Figure [51 where we show the error 
and constraint energies for r = 7 for the three evolution 
systems considered here. 

The results with the ADM formulation are drastically 
different than the corresponding results in the ID case. 
In the gauge condition we employed and in ID the ADM 
evolution equations can be shown to be strongly hyper- 
bolic, but in 2D they are only weakly hyperbolic. There- 
fore, due to the ill-posed nature of the ADM formulation 
the simulations of the 2D gauge wave crash quickly, and 
most importantly they crash faster with increasing reso- 
lution. We find that for r = 5, 6, 7 the ADM runs crash 
after 16,9,5 light crossing times respectively. 

The KST evolution equations are strongly hyperbolic, 
and hence the KST runs of the 2D gauge wave are longer 
than those of the ADM system. However, the simula- 
tions with the KST system crash very early on and can- 
not complete more than 30 crossing times. We have not 
been able to understand the mechanism via which the 
KST formulation crashes in this test. Similarly to the 
ID case, there are C^ij constraints which are excited be- 
cause of violation the momentum constraints, and then 
grow exponentially with time . However, when the com- 
putations stop the values of C^ij are of the order of 10~^° 
and this level of violation is not strong enough to explain 
the unexpected behavior of the KST system. The viola- 
tions of the momentum constraints introduce significant 
error into the Cuj and C2ij, but all constraints seem to 
grow exponentially at the same time, and thus we have 
not been able to attribute the blow up of the simulations 
to same source as in the ID gauge wave. Finally, we note 
that the instability which causes the evolutions to crash 
occurs at a time-scale which is resolution- independent. 

The results of the gauge wave simulations with the 
PADM formulation are drastically different. Figures [5] 
show that PADM successfully completes 100 light cross- 
ing times. In addition to that and as expected, the 
PADM formulation has better control of the constraints 



than the ADM and KST formulations, see Figure 9(a) 



Figure 9(b) shows that the PADM error energy experi- 
ences rapid growth. The same is true for the ADM and 
PADM error energies until the ADM and KST runs crash. 
This growth results because of excitation of exponential 
harmonic gauge modes as in the ID gauge wave. 

Finally, we note that for resolutions r = 5 and r — 6 
the PADM system manages to evolve the 2D strong gauge 
wave initial data for 56 and 98 light crossing times respec- 
tively. We have not been able to understand the mech- 
anism via which the simulations of the 2D gauge wave 
with the PADM system crash for low resolutions, but it 
seems to be correlated with the exponential growth of 
the solution. Nevertheless, the indisputable conclusion 
of the gauge wave simulations is that the numerical per- 
formance of the PADM system is much better than those 
of the KST and ADM formulations. 
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FIG. 9: Errors as functions of coordinate time t in units of light crossing times for the strong two-dimensional gauge wave and 
resolution r=7. On the left panel (a) we show the constraint energies and on the right panel (b) we show the error energies for 
the ADM, KST and PADM formulations. 



C. The linear wave test 

The linear wave test is a traveling gravitational wave 
of amplitude small enough so the evolution remains in 
the linear regime. The metric is given by 

ds'^ = -dt^ + dx^ + (1 + b)dif + (1 - b)dz^, (59) 

where 

fe = ylsinf— ^-^ — ^j, (60) 

where A is the amplitude of the wave and d the size of the 
evolution domain. The suggested amplitude for this test 
is ^ = 10"^ [10]. For the first-order ADM formulation 
((59| , deO]) is equivalent to 

711 = Ij 722 = 1 + fe, 733 = 1 - 6, (61) 




and all other dynamical variables are zero. 

Equation ([59)1 suggests geodesic slicing for the lapse 
function, but the test can be run both with a. — 7^/^ and 
a = 1 + ln7. This is so because for (|59p 7 = 1 — 6^. 
However, 6 is so small that a = 7^/^ « 1 — 6^/2 w 1 — 5 • 
10~^^ « 1. Thus, because of roundoff error a computer 
cannot "see" the difference between geodesic slicing and 
harmonic slicing. Similarly, it can be shown that "l-|-log" 
slicing is compatible with (|59p . As in the gauge wave 



case, the 2D linear wave can be obtained from ([5^ by 
the coordinate transformation fM]) . 

We carry out simulations of both a one-dimensional 
and a two-dimensional wave with harmonic slicing a — 
7^/^ and resolutions r = 5,6,7. The simulation time is 
1000 light crossings in the ID case and 100 light crossings 
in the 2D case. The results of our runs are summarized in 
Figures [TUl and [TT] where we show the constraint energies 
for the three formulations in the ID and 2D version of the 
test respectively. For brevity we do not show the error 
energies, because all formulations can easily handle this 
test and the essential difference between them lies in the 
behavior of the constraint energies. 

In the ID simulations (Figure [TU)) the constraint en- 
ergies of the KST and ADM formulations overlap and 
stay constant throughout the 1000 light crossings of the 
simulation. On the other hand, the PADM formulation 
pushes the constraint energy down close to the roundoff 
error level. There is a difference of at least 3 orders of 
magnitude between the constraint energy of the PADM 
system and those of the KST and ADM systems. 

In Figure [TT] we show the constraint energies of the 
three formulations considered in this work for the 2D lin- 
ear wave test. In the case of the PADM formulation we 
output the data at every light crossing time. However, 
outputting the constraint energies for the KST and ADM 
formulations every light crossing time, results in a false 
picture of the behavior of the KST and ADM formula- 
tions for this test, as the constraint energies can be seen 
to grow in time. What really happens is that the con- 
straint energies wildly oscillate roughly around the value 
2.4 • 10^^° and they show no sign of growth in time. For 
this reason, instead of outputting the constraint energy 
at every light crossing time, we output at every time- 
step and average the constraint energy through each light 
crossing. Then this average value is assigned to be the 
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FIG. 10: Constraint energies as functions of coordinate time 
t for the ID linear wave for resolution r = G. 
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FIG. 11: Constraint energies as functions of coordinate time 
t for the 2D linear wave test for resolution r = 6. 



value of the constraint energy at every light crossing time, 
which is what we have plotted in Figure [Til for the KST 
and ADM systems. 

From this figure it is clear that the ADM and KST 
constraint energies stay roughly constant throughout the 
entire run. For the 2D high amplitude gauge wave the 
damping properties of the PADM formulation are still 
evident, but not at the same level as in the ID case. 
Nevertheless, the PADM system controls the constraint 
violations better than the KST and ADM systems. 



D. Gowdy Spacetimes 

The polarized Gowdy T3 spacetimes are solutions of 
the Einstein equations which describe an expanding (or 
contracting) universe containing plane polarized gravita- 
tional waves, see for example [10.] and references therein. 



1. Expanding solution 
The expanding Gowdy metric is usually written as 



ds^ = t-^''^e^''^{-dt^ + dz^) +tdw^ , (64) 



where 



dw"^ = e^dx^ + e ^ dy^ . 



(65) 



The solution which is considered standard is the one 
with 

A(z,t) =27r2i2 (Jo(27rt)2 + Ji{2Titf) 

- 27rtJo(27ri)Ji(27rt) cos^(27rz) 

+ 7rJo(27r) Ji(27r) - 271^ (Jo(27r)2 + Ji(27r)2) 

(66) 



and 



P(z, t) = Jo(27rt) cos(27r.z), 



(67) 



where J„ denotes the Bessel function of the first kind of 
order n. 

The non-zero dynamical variables in the context of the 
first-order ADM formulation can be directly derived from 
([64)1 and they are presented in appendix [X] By use of 
equations ([64|) and ([65]) . it is straightforward to show 
that the lapse is harmonic and that 



(68) 



The expanding Gowdy wave test has proved to be one 
of the most challenging problems in numerical relativity. 
This is because the metric components grow exponen- 
tially with time and hence the truncation errors intro- 
duced in the numerical solution grow with time, if the 
resolution is fixed. Consequently, the numerical solution 
soon veers off the analytical one. Another result of the 
exponential growth is that very soon the dynamical vari- 
ables grow so large that the computations cannot be han- 
dled by using standard 64-bit floating-point arithmetic. 
Therefore, all simulations of the expanding Gowdy solu- 
tion are expected to be terminated [l3|. 

We start our expanding gowdy wave simulations ait = 
1 and evolve the initial data forward in time using the 
densitized lapse We run the test for resolutions 

r — 6, 7, 8 and until the code crashes. 

The results of the simulations are summarized in Fig- 
uresUll where we plot t he cons traint energy (Fig. 12(a) ) 
[T2(b)| ) for the ADM, KST and 
ution r = 8. These plots show 



and error energy (Fig 
PADM systems for reso 



that the runs with the KST formulation explode very 
quickly, whereas the ADM and PADM formulations can 
ran the test for significantly longer times. The error en- 
ergy and constraint energy of the PADM system lies lower 
than that of the ADM system and for this reason the 
PADM formulation extends the lifetime of the simula- 
tion by 12 light crossing times. 
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FIG. 12: Errors as functions of coordinate time t for the expanding Gowdy wave for resolution r = 8. On the left panel (a) we 
show the normalized constraint energies and on the right panel (b) we show the normalized error energies for the ADM, KST 
and PADM formulations. 



We checked that the simulations with the KST sys- 
tem do not crash because of violation of the numerical 
stability criterion. Instead, like in the gauge wave and 
robust tests, the error growth of the Gowdy simulations 
with KST system can be explained by the mixing and 
excitation of unwanted degrees of freedom the KST equa- 
tions cause, because of violations of the momentum con- 
straints. However, the situation with the expanding po- 
larized Gowdy wave is much more severe than the gauge 
wave. This is simple to understand by a glimpse at the 
evolution equation of C333 which is 



dtC: 



'tl--333 



r 0733-^3, 



(69) 



Equation ((69|) shows that, because of the exponential 
growth of both the lapse function and the three-metric 
with time, the C333 constraint is bound to grow expo- 
nentially with time. The RHS of ([55]) varies roughly as 
^-3gA/4^^^ By using the analytic solution we can esti- 
mate the RHS of For example, at t = 50 even if 
Ms were M3 = 10"^^ the RHS of equation ^ is of the 
order of 10^'^. 

Although we are able to evolve expanding Gowdy wave 
initial data with the ADM and PADM systems to very 
late times, it is seen in Figure 12(b) that at the highest 
resolution the numerical solution matches the analytical 
solution up to approximately 110 light crossings with the 
ADM formulation and up to 122 light crossings with the 
PADM formulation. At these times the normalized er- 
ror energy becomes of order 0.1. The accuracy with the 
PADM system is completely lost after roughly t = 275 
and with the ADM system after roughly t — 258. Wc 
note here that in [s^ the ID (not 3D) simulations of the 
expanding gowdy wave with the KST system completely 
lose accuracy after 150 light crossing times, even though 
the authors use a much more sophisticated and accurate 



method of integration (i.e., pseudo spectral methods) in 
order to perform their simulations. For resolution r = 9 
we can accurately evolve the solution with the PADM 
system up to approximately 200 hght crossings. 

Taking all these facts into accounts we conclude that 
the comparison of the three formulations in the expand- 
ing Gowdy wave testbed shows that the PADM formu- 
lation performs better than the ADM formulation and 
significantly better than the KST formulation. 



2. Collapsing Gowdy Cosmology 

Unlike the expanding solution, which moves away from 
the t — singularity, the collapsing solution will reach 
the singularity very soon if integrated in the coordinate 
time t as is dictated by the metric (|64p. For this reason a 
coordinate transformation is implemented by asking for 
a new coordinate time r, such that t = F(t) = ke'^'^ 
with k,c constants given by c ~ 0.002119511921460752 
and k ~ 9.670769812764059. In this new coordinate time 
the singularity is reached at r = — cxi and the numerical 
integration is much more tractable. The metric in the 
new coordinates is given by 



where 



and 



i(r) = cF{t) 



3/4gA(^,F(T))/4 



(70) 



(71) 



dw^ = e^(^'^(^»dx2 -f e-^(-^^M)dy2. (72) 

The non-zero dynamical variables in the context of the 
first-order ADM formulation can be easily derived by use 
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FIG. 13: Errors as functions of coordinate time t for the collapsing Gowdy case for resolution r—7. On the left panel (a) we 
show the normalized constraint energies and on the right panel (b) we show the normalized error energies for the ADM, KST 
and PADM formulations. 



of ((7D|) and they are presented in appendix [Bl The lapse 
function is harmonic and it is straightforward to show 
that 



(73) 



We start the simulations at the 20th zero of the Jo{2TTt) 
Bessel function, i.e., at Tq ~ 9.87532058290983 and 
evolve the initial data backwards in time using ([75)1 . We 
run the test for r — 6, 7, 8 and for 1000 light crossing 
times with all formulations considered in this work. Fi- 
nally, we note that for this test the values of the PADM 
parameters are opposite to those of equation P7)) . 

The results of our simulations are summarized in Fig- 
ures [13] where we show the constraint energies and the 
error energies for all formulations considered here. We 
can see in 13(a) that the ADM constraint energy grows 



in time and becomes larger than the PADM and KST 
constraint energies very early on. At 1000 light crossing 
times the ADM constraint energy is about two orders of 
magnitude larger than the PADM and KST constraint 
energies. However, we see in Figure |13(b)| that the ADM 
error energy lies below the PADM and KST error en- 
ergies for up to about 300 light crossing times. After 
500 crossing times the ADM error grows faster than the 
PADM and KST error and at 1000 crossing times there 
is roughly one order of magnitude difference between the 
ADM error energy and the corresponding error energies 
of the PADM and KST formulations. 

Figures fT3l show that the PADM and KST systems per- 
form equally well. However, for the most part of the 1000 
crossing times the PADM constraint and error energies 
are less than those of the runs with the KST system. To- 
wards the end of the simulation the KST error energy 
is slightly less than the PADM error energy, while the 
PADM constraint energy is slightly less than the KST 
constraint energy. 



V. CONVERGENCE 

In this section we test the convergence and accuracy 
of our numerical schemes. In order to test the order of 
convergence we need to use exact solutions to the sets of 
PDEs which we solve. Of the four tests we performed for 
the comparison of the ADM, KST and PADM formula- 
tions, only the gauge wave and polarized Gowdy space- 
times are exact solutions of the Einstein equations and 
only these can form a real basis for convergence testing. 

Unfortunately, the gauge wave results suffer by the ex- 
citation of exponential harmonic gauge modes and for 
this reason they cannot serve as a basis for convergence 
testing for long time. This is so, because of the prop- 
erties of the solutions and not because of our numerical 
scheme. We demonstrated in section |IVB] that for differ- 
ent resolutions different exponential modes are excited 
and at different time-scales. Thus, convergence testing 
with the gauge wave can be done only for short periods 
of time before exponential modes become noticeable. We 
have been able to show that our code is indeed conver- 
gent and second-order accurate for short periods of time 
in the gauge wave test. 

The expanding Gowdy test cannot serve as a basis for 
convergence testing for long times either because the test 
cannot be run for very long times. For example the KST 
formulation cannot even complete 22 light crossing times 
in the expanding Gowdy cosmology. Again this is be- 
cause of the properties of the solution and not because of 
our numerical scheme. Nevertheless, we have been able 
to show that our code is convergent and second-order ac- 
curate in simulations of the expanding Gowdy wave with 
the ADM, KST, and PADM systems up until the sim- 
ulations crash. In order to keep this discussion short, 
we omit the convergence plots for the expanding Gowdy 
and gauge wave tests, but they can be available by the 



17 



-1 



-3 



- 




















- 
























'■^^ 

r=7,,- 






/■••/"'• 








it 


'.I'AV'f ' 


>..-•/'»"■'■ 


r=8,- 












- 


yfp ■ 









































-1 -2 -3 -4 -5 -6 -7 -8 -9 -10 
t/100 

FIG. 14; Normalized error energies as functions of coordinate 
time t, for resolutions r — 6,7, 8. The plot corresponds to 
the simulations of the collapsing Gowdy spacetime with the 
ADM formulation. 
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FIG. 15: Normalized error energies as functions of coordinate 
time t, for resolutions r = 6,7, 8. The plot corresponds to the 
simulations of the collapsing Gowdy spacetime with the KST 
formulation. 





ADM 


KST 


PADM 


5,6 


1.9057 


2.0461 


1.9617 


6,7 


1.9793 


2.0146 


1, 9945 


7,8 


2.0296 


2.0032 


2.0144 



TABLE I: Convergence order table. The first column shows 
the resolutions which are involved in the calculation of the 
order of convergence and the second to fourth columns show 
the order of convergence of our code with the ADM, KST and 
PADM formulations respectively. 

authors upon request. 

In order to test long term convergence of our code, we 
will use our results of the collapsing Gowdy wave testbed, 
where all three formulations successfully complete 1000 
light crossing times. The convergence plots can be seen 
in Figures [H [15] and [16] for the ADM, KST and PADM 
formulations respectively, where it is evident that with 
increasing resolution the numerical solution converges to 
the analytical one for all formulations. 



It is also important to demonstrate that the order of 
convergence of our code is 2. To do this we use our data 
for the error energies and calculate the order of conver- 
gence which is given by 

T{tn,ri,r2) = log2 II „, 7|T ' 

r2-ri \\\5UN{tn,r2)\\J 

where ri and r2 two arbitrary resolutions, tn a given 
time-step and ||(5Wjv(tn, r-)!! is the normalized error en- 
ergy. Exact second-order accuracy means that T ~ 2. 
In table H] we show the time-average of T for the col- 
lapsing Gowdy spacetime with all three formulations 
and resolutions {ri,r2} = {5,6}, {ri,r2} = {6,7} and 



{ti, ''2} = {7, 8}. The time-average order of convergence, 
T, is 

- 1 ^ 

^{ri,r2) = —^T{tn,ri,r2), (75) 

n=Q 

where in our case tn — n and N — 1001. Tableland 
Figures [TIKTB] clearly demonstrate that our code is con- 
vergent and second-order accurate. 
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FIG. 16; Normalized error energies as functions of coordinate 
time t, for resolutions r = 6, 7, 8. The plot corresponds to 
the simulations of the collapsing Gowdy spacetime with the 
PADM formulation. 



VI. CONCLUSIONS 

We have described a stable, convergent and second- 
order accurate numerical scheme for solving the equa- 
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tions of the recently proposed PADM formulation of GR, 
we have tested the accuracy and stability of the PADM 
system and compared it with the first-order ADM and 
KST evolution systems. 

The ADM system is generally ill-posed, but in con- 
jmiction with a densitized lapse or "1+log" slicing the 
ID ADM evolution equations are strongly hyperbolic and 
can be compared to other well-posed systems. The KST 
evolution system in conjunction with either a densitized 
lapse or "1+log" slicing is strongly hyperbolic. The well- 
posed PADM equations have the structure of a mixed 
hyperbolic second-order parabolic set of PDEs. PADM 
possesses the desirable property of constraint damping. 
Unlike all other systems used in numerical relativity with 
constraint damping, the PADM system damps the short 
wavelength constraint violating modes more efficiently 
than the long wavelength ones. 

The scries of tests which we used to compare the ADM, 
KST and PADM systems is known as the "Apples with 
Apples" or the "Mexico City" tests. These tests are 
designed to probe both the linear and the non-linear 
regime of the Einstein equations and they are called the 
robust stability test, the gauge wave stability test, the 
linear wave stability test and the Gowdy wave stability 
test. The first three can be considered as perturbations 
about flat spacetime, so they probe the weak field limit 
of the Einstein equations, whereas the last one probes 
the strong field regime. 

We ran the robust stability test with a densitized lapse 
with a = 0.5, 1 and "l-|-log" slicing and we demonstrated 
that the PADM formulation performs much better than 
both the ADM and KST systems. The constraint vi- 
olations with PADM are about seven orders of magni- 
tude smaller than the constraint violations with KST 
system. PADM damps the constraints very quickly and 
the shorter the wavelength of the constraint violations 
the faster they are damped. This is in complete agree- 
ment with the theoretical properties of PADM. Finally, 
the ADM formulation crashes very quickly and faster for 
higher resolutions. Finally we demonstrated that the be- 
havior of all three systems in this test is similar in the 
other gauge conditions we employed. This confirms the 
theoretical prediction that the PADM formulation damps 
the constraints independently of the gauge chosen. 

The simulations of the high amplitude ID gauge wave 
test clearly demonstrate that the PADM formulation per- 
forms better than the ADM formulation, while the KST 
system cannot evolve the initial data for more than 200 
crossing times. The reason for this behavior of the KST 
system is that the dynamics of the formulation is such 
that it excites imwanted degrees of freedom, due to error 
in the momentum constraints. These degrees of freedom 
are not part of the dynamics of the gauge wave testbed 
and after they are excited, they grow exponentially with 
time and lead the computations to an end. 

The simulations of the high amplitude 2D gauge wave 
test also demonstrate that the numerical performance of 
the PADM formulation is much better than those of the 



ADM and KST systems. In this case it is the ADA/[ simu- 
lations which crash very quickly. The runs with the KST 
equations cannot evolve the initial data for more than 
30 crossing times no matter how large the resolution is. 
Finally, the runs with PADM can successfully complete 
100 light crossing times and show that the PADM system 
damps the constraints with time. 

The linear wave testbed presented no essential diffi- 
culty for any of the three formulations, but both the ID 
and the 2D linear wave simulations, show that the PADM 
formulation performs better than the ADM and KST for- 
mulations. In the ID case the constraint error with the 
PADM system is quickly damped and is three orders of 
magnitude less than the constraint error with the ADM 
and KST systems. We observed similar behavior in the 
2D linear wave, but the difference in the constraint error 
among the three formulations is not as large. 

The polarized Gowdy wave testbed has two versions: 
the expanding one and the collapsing one. In the collaps- 
ing Gowdy wave, the KST and PADM systems perform 
equally well, but much better than the ADM one. At the 
end of the simulations the ADM constraint error is at 
least two orders of magnitude larger than the constraint 
error of the PADM and KST systems and the error in the 
solution with the ADM system is one order of magnitude 
larger than the corresponding errors with the KST and 
PADM systems. 

The expanding Gowdy wave was by far the strongest 
field solution we simulated. The crash of simulations of 
the expanding Gowdy solution, is naturally expected for 
any evolution system because the metric grows exponen- 
tially with time. The results of this test show that the 
PADM formulation performs better than the ADM for- 
mulation and dramatically better than the KST formula- 
tion. The expanding Gowdy wave initial data cannot be 
evolved with the KST system for more than 25 crossing 
times at the highest resolution we attempted. The cause 
for the sudden break down of the KST system is the same 
as the cause for the sudden break down of the KST sys- 
tem in the ID gauge wave. In contrast, the ADM and 
PADM systems can evolve the expanding Gowdy wave 
solution for significantly longer times, until the error in 
the solution grows so large that simulations terminate. 
The high resolution ADM runs completely lose accuracy 
after 258 light crossing times. The PADM runs are more 
accurate than the ADM ones and at the highest reso- 
lution completely lose accuracy after 275 light crossing 
times. 

Taking all these facts into accounts, we conclude that 
PADM successfully passes the standard tests of numer- 
ical relativity and works equally well with a variety of 
algebraic gauges. Via the comparison of the numeri- 
cal performance of the PADM formulation and those of 
the ADM and the KST formulations, we conclude that 
PADM has better control of the constraint violations 
than both ADM and KST. The PADM system performs 
better than the ADM one and in most tests better than 
the KST system. 
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APPENDIX A: EXPANDING GOWDY 
SOLUTION 



directly derived by use of equations ([701) -([111), ^ and ^ 
and they are given by 



In this appendix we present the non-zero dynamical 
variables, in the context of the first-order ADM formu- 
lation, for the expanding Gowdy solution. Those can be 
directly derived by use of equations ([64ll - (|67|) . Q and ([2]) 
and they are given by 



te~^, 7,,=t-V2eA/2^ (Al) 



and 



Iz. = F{t) 



-l/2gA(^,_F(T))/2 



(Bl) 



where 



K 



yy 



2a 
1 

2^ 



4a 



a = V72 



(A2) 



(A3) 



c dP 

Kyy--^lyya^F^), 



(B2) 



where the lapse function a is given by (|7ip . In addition, 
we have 



dtX{z,t) = 27r2t[(cos(47rz) + l)Ji{2Trt)^ 

+ 2Jo(27rt)^sin2(27rz)] 



and 



(A4) 



dtP{z,t) = -27rJi(27ri)cos(27rz). (A5) 
We also have 

(A6) 
(A7) 



where 



dzP{z, t)^- 27rJo(27ri) sin(27rz), 
dz\{z,t) = 47r2iJo(27ri)Ji(27ri)sin(47rz) 



^ 27r2F[(cos(47rz) + l).U2^Ff ^^^^ 
+ 2jQ{2T:Ffsm^{2 



and 



dP{z, F) 



-27rJi(27ri^)cos(27rz) 



(B4) 



Finally the non-zero components of the Dkij variables 
are given by 



DzXX = IxxdzP, Dzyy = ^-fyyOzP, D zzz = ■7pZzdz\ 



(B5) 



where 



APPENDIX B: COLLAPSING GOWDY 
SOLUTION 

In this appendix we present the non-zero dynamical 
variables, in the context of the first-order ADM formu- 
lation, for the collapsing Gowdy solution. Those can be 



dzP{,z, F) = - 27tJq{2ttF) sin(27rz), 
dz\{z,F) = 47r2FJo(27rF)Ji(27rF)sin(47rz) 



(B6) 
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